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I. INTRODUCTION 

The aim of this paper is to present a computation of the 
interfacial tension 71N between the coexisting isotropic 
and nematic (IN) phase in suspensions of monodisperse 
hard rods via computer simulation. While the hard-rod 
fluid simplifies experimental reality, ignoring for example 
long-ranged interactions and polydispersity [l|, it nev- 
ertheless captures the main mechanism of the IN phase 
transition and serves as a valuable model system. Exper- 
iments have shown that 7™ is very small, typically in the 
range 10 _3 -10" 4 mN/m ,1], which makes it difficult to 
extract 71N from simulation data. Simulation estimates 
of 71N are therefore rare, and have only been reported for 
ellipsoids 0, Q , soft rods 3 , and lattice models[3L The- 
oretical estimates are more abundant [1 S 11 |SM El , 
but are usually obtained in the Onsager limit 12] of in- 
finite rod length (L/D — > 00). The case of finite rod 
length is more difficult to describe theoretically, but has 
been addressed in Ref. |9| using density functional theory, 
and in Ref. H using a scaling relation. At the time of 
writing, no simulation estimate of 71N for the hard-rod 
fluid has been reported. Such an estimate would clearly 
be valuable to test theoretical predictions, and to see if 
the order of magnitude of 71N observed in experiments is 
reproduced. 

Despite its simplicity, simulating the hard-rod fluid is 
not trivial 0, El . The bottleneck is the hard particle 
interaction, which complicates both molecular dynamics 
(MD) and Monte Carlo (MC) methods. In the case of 
MD, the discontinuous potential prevents the calculation 
of smooth forces. In the case of MC, equilibration times 
are long due to very low acceptance rates. An important 
improvement is the use of soft interactions, as was done 
for ellipsoids 0,0 and rods 0,E1- By using soft interac- 
tions, the qualitative phase behavior is usually retained, 
but simulations become much more efficient. Moreover, 
MC simulations in the grand canonical ensemble become 
possible, enabling the investigation of IN coexistence via 
the probability distribution in the particle number den- 
sity. This technique is well established in simulations of 
fluid- vapor coexistence 0, El fl8L IT^. I20I ] , and was re- 



cently extended to IN coexistence in suspensions of soft 
rods 0. The advantage of grand canonical simulations 
is that the coexistence densities, as well as the interfacial 
tension, can be obtained. 

Since coexisting phases are separated by a free energy 
barrier arising from the interfacial tension |2lj . it is es- 
sential to use a biased sampling scheme to access regions 
of high free energy. In simulations of fluid-vapor coex- 
istence, the bias is usually put on the density. While a 
density bias has also been used to simulate IN coexis- 
tence 0, this choice is not optimal. In simulations that 
rely on standard MC moves, such as random translations 
and rotations of single particles, it is difficult to reach the 
nematic phase starting in the isotropic phase simply by 
increasing the density because the orientational degrees 
of freedom relax only very slowly [2^ • This effect is called 
"jamming" , and it explains why the simulations of Ref. 4 
were limited to rather small systems. 

In this work, grand canonical MC simulations using a 
bias on the nematic order parameter are performed. As 
we will show, this approach is much less susceptible to 
jamming, and enables simulations of large systems. This 
in turn allows for accurate estimates of the interfacial ten- 
sion in suspensions of soft rods. As an additional bonus, 
a bias on the nematic order parameter paves the way to- 
wards grand canonical simulations of hard rods, enabling 
a simulation estimate of 71N for the hard-rod fluid. 

The outline of this paper is as follows: First, we in- 
troduce the liquid crystal model used in this work. The 
biased sampling scheme is described next. The results 
are presented in Sec. IIVI We end with a summary and 
comparison to theoretical predictions in the last section. 

II. MODEL AND ORDER PARAMETERS 

We consider rods of elongation L and diameter D. 
The simulations are performed in a three dimensional 
box of size L x x L y x L z using periodic boundary con- 
ditions in all dimensions. In this work, we fix L x = L y , 
but we allow for elongation in the remaining dimension 
L z > L x . Moreover, to avoid double interactions between 
rods through the periodic boundaries, we set L x > 2L. 
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The position of the center of mass of rod i is denoted f^, 
and its orientation Wj, with normalization |u,-| = 1. The 
interaction between two rods i and j is given by a pair 
potential of the form 



(r) 



e r < D, 
otherwise, 
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with r the distance between two line segments of length 
L, see also Ref. 0. The total energy is thus a function of 
the center of mass coordinates and the orientations of all 
rods 
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with N the number of rods in the system (in the following 
we will drop the fi and Ui dependences in our notation). 

To investigate the IN transition, the density and the 
average rod alignment are used as order parameters. 
Since the density in the isotropic phase is lower than 
in the nematic phase, the rod number density p = N/V 
may be used to distinguish the phases, with V the vol- 
ume of the simulation box. Following convention, we 
also introduce the reduced density p* = p/p cp , with 
p cp = 2/[v / 2 + (L/D)^/S\ the density of regular close 
packing of hard rods [T]| . 

In the nematic phase the rods are on average aligned, 
whereas in the isotropic phase the rods are randomly ori- 
ented. Therefore, the nematic order parameter may also 
be used to distinguish the phases. The latter quantity 
is defined in terms of the orientational tensor Q, whose 
components are given by 
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with Ui a the a component (a = x, y, z) of the orientation 
of rod i, and 5 a p the Kronecker delta. In this work, 
the maximum eigenvalue S of the orientational tensor is 
taken as nematic order parameter, being close to unity 
in the nematic phase, and close to zero in the isotropic 
phase. The eigenvector corresponding to S is called the 
director, and it measures the preferred direction of the 
rods in the nematic phase. 



III. SIMULATION METHOD 

We study IN coexistence via grand canonical MC sim- 
ulations. In the grand canonical ensemble, the volume, 
the temperature T, and the chemical potential p are 
fixed, while the number of rods in the simulation box 
fluctuates. Insertion and removal of rods are attempted 
with equal probability and accepted with appropriate 
Metropolis rules to be given later. The aim of grand 
canonical simulations is to measure the probability dis- 
tribution in the number of particles P(N). At the coex- 
istence chemical potential, P(N) becomes bimodal with 



two peaks of equal area. An example distribution is 
shown in Fig. 0] where we have plotted the logarithm 
of P(N). The peak locations yield the coexistence densi- 
ties; the average height of the peaks Af2 in fcBTlnP(TV) 
is the free energy barrier separating the phases, with fee 
the Boltzmann constant. In three dimensions using pe- 
riodic boundary conditions and for sufficiently large sys- 
tems, the barrier is related to the interfacial tension via 
7in = AJ7/(2L^), where L x is the lateral dimension of 
the simulation box [23| . 

In simulations, the free energy barrier presents a prob- 
lem. Unless Aft is small, such as close to a critical point, 
simulations rarely cross the barrier, and spend most time 
in only one of the two phases. Biased sampling tech- 
niques are required to overcome the barrier. In general, 
these techniques aim to construct a weight function W(Q 
of some bias variable C- The weight function is con- 
structed such that a simulation using a modified potential 
E'(() = E + fceT VF(C) yields a uniform probability dis- 
tribution in the bias variable, with E the potential of the 
original system. The grand canonical acceptance rules 
using the modified potential read as 
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for the insertion and removal of a single particle, respec- 
tively j3)E3- Here, Co and Ci denote, respectively, the 
value of the bias parameter in the initial and final state, 
AE is the potential energy difference between initial and 
final state given by Eq.J2J), and /3 = l/{k^T). For a 
properly constructed W(£), the biased simulation sam- 
ples all states ( with uniform probability. Once W{C) 
is known, the distribution P(N) of the unbiased system 
can be constructed. 

One is rather free in choosing the bias variable. The 
best choices are variables that change significantly when 
going from one phase to the other. For fluid- vapor transi- 
tions, a natural bias is the particle number density. In the 
case of IN coexistence, the density is still a valid variable 
because of the density gap between the isotropic and the 
nematic phase. This was used in Ref. to study IN co- 
existence in suspensions of soft rods. Whether a bias on 
density in systems of elongated particles is efficient, de- 
pends on how easily a dense isotropic phase can rearrange 
itself to become nematic. In practice, the jamming effect 
limits density biased sampling to rather small systems 
and soft interactions. As it turns out, for IN transitions, 
a much more powerful bias variable is the nematic or- 
der parameter 5*. Note, however, that phase coexistence 
is defined in terms of P(N). Therefore, in a simulation 
which biases on S, the distribution P(N) must still be 
reconstructed. To this end, histograms in both the par- 
ticle number N, as well as in S, have to be measured. In 
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this section, we explain how the bias on S, and the sub- 
sequent reweighting in N and S, are implemented. It is 
convenient to discuss the more straightforward procedure 
of a density bias first. 

A. Biased sampling on p 

A convenient method to bias on density is Successive 
Umbrella Sampling (SUS) |26| . Here we describe the al- 
gorithm in its simplest form; refinements are given in 
the original reference. The choice of the sampling algo- 
rithm is not crucial. The general principles also apply to 
other schemes, such as conventional umbrella sampling 
|27j . multicanonical sampling j2^|, Wang -Landau sam- 
pling 29] , or hyperparallcl tempering 30] . 

The aim is to construct a function W(N) of the num- 
ber of particles such that a simulation using the mod- 
ified potential E'(N) = E + k B TW(N) yields a uni- 
form distribution in N, with E given by Eq.J2J). The 
modified potential thus contains an explicit dependence 
on the bias variable N . Following Ref. I26L the particle 
number axis is divided into equally sized intervals called 
windows, starting with some minimum number of par- 
ticles Nq. In the first window, the number of particles 
is confined to Nq < N < iVo + 1, in the second window 
to Nq + 1 < N < No + 2, and in the i-th window to 
Nq + i — 1 < N < Nq + i. In this example, the win- 
dow size equals a single particle but this choice is not 
essential: SUS works just as well using larger windows 
|26| . The choice of the window size is not completely ar- 
bitrary. Choosing the windows too large leads to poor 
sampling statistics at the window boundaries; choosing 
the windows too small runs the risk that certain relax- 
ation pathways are cut-off. In practice, a compromise 
needs to be made. 

The idea of SUS is to construct W(N) by simulating 
the windows separately and successively. Starting in the 
first window (i — 1), grand canonical MC moves are at- 
tempted (optionally combined with canonical moves such 
as translations and rotations), with the constraint that 
states outside the window bounds are rejected to fulfill 
detailed balance at the window boundaries. The relevant 
weights in the first window are W(N ) and W(N + 1), 
which we initially set to zero. We then record fl and 
counting the occurrence of the state with Nq and Nq + 1 
particles, respectively. In this notation, the subscripts 
"L" and "H" refer to the "lower" and "higher" window 
bound, respectively, while the superscript refers to the 
window number. To obtain a uniform distribution in N, 
the ratio of the counts should be unity. This will gen- 
erally not be the case, but is enforced by updating the 
weight of the higher window bound to 

W ncw (No +i) = W old (No + i) + Hfli/ft), (6) 

leaving the weight of the lower bound W(Nq + i — 1) 
unchanged, where i is the window number. In case > 
/£, the effect of this modification is a lower insertion 
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FIG. 1: Schematic illustration of biased sampling on the 
nematic order parameter. See details in text. 



rate, see Eq.Q, and a higher removal rate, see Eq.JSJ, 
leading to a count ratio closer to one. The latter can 
be checked by performing a second simulation using the 
updated weight (the reasoning for < /£ is similar). In 
practice, it may occur that one of the counts is zero. It is 
then necessary to modify W(No + i) by hand first, before 
starting the simulation. Note also that long simulation 
runs may be required to obtain the count ratio accurately. 

Having simulated the first window, W(Nq) and 
W(N + 1) are known. The choice W(N ) = is ar- 
bitrary but has no physical consequences since it merely 
shifts the potential by a constant. Next, we consider 
window 2, where the number of particles is allowed to 
fluctuate between N + 1 and Nq + 2, with respective 
weights W(N + 1) and W(N + 2). An important opti- 
mization of Ref. 26> is to linearly extrapolate the known 
weights W(iVo) and W(N + 1) to obtain an estimate 
for W(iVo + 2) (note that for the third and subsequent 
windows, quadratic extrapolation can be used). The sim- 
ulations in the second window are then performed using 
the extrapolated estimate, and the respective counts, f£ 
and of visiting the state with Nq + I and Nq + 2 
particles, are recorded. Finally, the weight W(Nq + 2) is 
updated using Eq.JSJ, leaving the other weight W(Nq + 1) 
unchanged, and the next window is considered. 

The above procedure is repeated until all windows 
of interest have been simulated, and the corresponding 
weight function W(N) is constructed. The sought-for 
distribution in the number of particles P(N) is trivially 
obtained via P(N) = Ce w{N \ with C a normalization 
constant pi) . 



B. Biased sampling on S 

Next, we consider the extension to a bias on the ne- 
matic order parameter S. Here, the modified potential 
reads as E'(S) = E + k B TW(S), with E given by Eq.©. 
The aim is to construct W(S) such that a simulation us- 
ing the modified potential samples all values of S with 
uniform probability. 
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The windows are obtained by dividing the nematic or- 
der parameter axis into equally sized intervals of width 
AS*. In the first window, the nematic order parame- 
ter is confined to < S < AS, in the second window 
to AS/2 < S < 3AS/2, and in the i-ih window to 
(i - l)A5/2 < S < (i + l)AS/2, see Fig. Q] The win- 
dows thus partially overlap. To sample both the isotropic 
and the nematic phase, the sampling range should span 
from S = to S « 1. Note that S is a continuous vari- 
able, whereas the density (expressed in the number of 
particles) is discrete. Therefore, a natural width for the 
windows does not exist, and one is forced to choose AS 
rather arbitrarily. We found that AS « 0.001 - 0.002 
gives good results, which means that O(10 3 ) windows 
are required to sample the transition. A consequence of 
discretizing the nematic order parameter is that W(S) is 
defined in steps of AS/2. Therefore, in the i-th window, 
W(S) assumes only two distinct values 

W(S) = i Wi ~ l — \)AS/2 < S < Sm 
[ ' [Wi S M <S< 0'+l)AS/2, 

with Sm = iAS/2 the center of the window (note that 
i>0). 

Starting in the first window (i = 1), the relevant 
weights are Wo and W\, which are initially set to zero. 
While simulating the first window, we count the occur- 
rence of states with < S < S M (ft) and S M < S < AS 
(/h)j with Sm = AS/2. To obtain the distribution in the 
number of particles P(N) (after all the quantity of inter- 
est) particle number histograms must also be stored (note 
that N fluctuates freely in each window). In the first win- 
dow, we thus record the probability distribution in the 
number of particles p^(N) for states with < S < Sm, 
and p]i(N) for states with 5m < S < AS. It is rec- 
ommended to store the distributions unnormalized. This 
makes it more convenient to restart the simulations at 
a later stage in case higher precision is required. After 
simulating the first window, the weight of the higher win- 
dow bound is updated to force a uniform distribution in 
S using W^ n cw = ^,oM + M/h//l), while keeping Wi-i 
fixed, where i is the window number. 

In the second window (i = 2), the relevant weights are 
W\ and W2- To simulate efficiently, the weights Wq and 
W\ of the previous window are extrapolated to estimate 
W2 ■ The extrapolated estimate is used while simulat- 
ing the second window, and the counts and are 
recorded, as well as the distributions pf,(N) and p^(N). 
After simulating the second window, Wi is updated as 
before, and the next window is considered. 

The above procedure is repeated up to some maximum 
number of windows w max chosen well into the nematic 
phase. The remaining step is to combine the weights 
Wi with the distributions p\,(N) and pfj(iV) to obtain 
P(N). Note that the upper region of window i over- 
laps with the lower region of the next window i + 1, 
see Fig. n More precisely, the distributions p^N) and 
p z l ^ 1 (N) stem from the same S interval and are thus 



measured with the same probability by the sampling 
scheme. Therefore, these distributions may be combined 
pi(N) = Pn(N) + p t l / 1 (N), and normalized such that 
Ew=oP'W = 1- The distribution P(N) is simply a 
weighted sum of the above (normalized) pi{N). Since 
— ksTWi corresponds to a free energy, each pi(N) con- 
tributes to P{N) with a weight proportional to e \ 
This leads to P(N) = C Yh^T Pz(N)e w % where the 
sum is over all windows, and normalization constant 

C. Bias on p versus bias on S 

Clearly, the discussed methods serve the same purpose: 
to measure the distribution P(N) at coexistence. Density 
biased sampling is by far the easiest to implement. It has 
the additional advantage that the coexistence chemical 
potential need not be specified beforehand: once P(N) 
has been measured at some chemical potential fio, it can 
be extrapolated to any other chemical potential \i\ by 
using the equation 

P(N\vi) = P(N\^ )e^-^ N , (7) 

with P(N\/i a ) the probability distribution P(N) at 
chemical potential fi a . Obviously, one should establish 
roughly beforehand the density at which the IN transi- 
tion occurs, to avoid sampling large regions of irrelevant 
phase-space. 

The situation is reversed when biasing on the nematic 
order parameter. In this case, the sampling range is al- 
ways from S = to S f=a 1. However, to observe phase 
coexistence, it is essential to use a chemical potential 
that is rather close to the coexistence value. Of course, 
Eq.(J7J) still holds, but the range in /i over which one can 
extrapolate is much smaller, precisely because the bias is 
put on S and not on p. An estimate of the coexistence 
chemical potential may be obtained in a density biased 
simulation of a small system, or via the Widom insertion 
algorithm [25l l3lj . This certainly makes biasing on S 
more involved. Moreover, for each attempted MC move, 
S in the final state must be determined, regardless of 
whether the move is accepted. It is therefore important 
to calculate S efficiently. In particular, the O(N) loop of 
Eq.© should be eliminated, which can be done following 
the method outlined in Ref. 



IV. RESULTS 

An important conclusion of Ref. 4 is that the IN in- 
terfacial tension obtained from P(N) may be prone to 
strong finite size effects. Away from any critical point, 
interfaces are the dominant source of finite size effects. 
The use of periodic boundary conditions leads to the for- 
mation of two interfaces. In small systems, the interfaces 
may interact and this will influence the estimate of 7in- 
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A convenient way to suppress interface interactions is to 
use an elongated simulation box with L z 3> L x |32| . in 
accord with Fig. This forces an orientation of the in- 
terfaces perpendicular to the elongated dimension (since 
this minimizes the interfacial area) , with a separation be- 
tween the interfaces that is larger than it would be in a 
cubic system of the same volume. The absence of in- 
terface interactions is manifested by a pronounced flat 
region between the peaks in lnP(iV). Note that a flat 
region is essential, but not sufficient, to extract 7in reli- 
ably. There may still be finite size effects in the lateral 
dimensions L x and L y , arising for instance from capillary 
waves. Ideally, the lateral dimensions should be large 
enough to capture the long wavelength limiting form of 
the capillary spectrum [8^ | . To actually measure the cap- 
illary spectrum of the IN interface is demanding [3J. A 
more convenient approach sufficient for our purposes is 
to first establish a minimum elongation L z in which in- 
terface interactions are suppressed, and to then check for 
finite size effects in the lateral dimensions by varying L x 
and L y explicitly. 

An additional motivation to use large lateral dimen- 
sions is to stabilize the interfaces. The interfacial free 
energy is of order jinL x , and if this is small compared to 
/cbT\ the interfaces will generally not be stable. These is- 
sues are especially relevant for IN coexistence because 71N 
is very small. Therefore, in this section, we first perform 
MC simulations in the canonical ensemble to obtain an 
indication of the system size required to observe stable 
interfaces. Next, we present coexistence data obtained 
using the nematic order biased sampling scheme. 



A. Interfacial profiles 

We consider hard rods, i.e. e — > 00 in Eq.Q, of elon- 
gation L/D = 15. The simulations are performed in 
the canonical ensemble, where the number of rods, the 
volume, and the temperature are fixed. The box dimen- 
sions are L x = L y — 10L/3 and L z = 20 L. We set 
the overall density of the system to p* = 0.205, which 
is well inside the coexistence region 0, ^| , correspond- 
ing to ca. 11, 000 particles. An initial system is prepared 
containing two interfaces, with the director of the ne- 
matic phase aligned in the plane of the interface. This 
is the stable configuration, as confirmed by theory 0, Q 
and simulation 0,0, El- The initial system is evolved 
with random rotations and translations of sing le rods, 
accepted with the standard Metropolis rules j2J, l25l l34j . 
The system is equilibrated for 10 6 sweeps, after which a 
snapshot is taken every 260 sweeps, up to a total of 3x 10 4 
snapshots (one sweep corresponds to one attempted MC 
move per rod). 

After equilibration, simulation snapshots schemati- 
cally resemble Fig. |5J Note, however, that they contain 
far more particles than depicted in this simple sketch. 
The aim is to measure the density profile p(z), and the ne- 
matic order parameter profile S(z) along the elongated z- 
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FIG. 2: Schematic representation of a simulation snapshot 
at IN coexistence. The ordered nematic phase is located in 
the middle of the simulation box. Profiles are measured along 
the ^-dimension using bin size Az. 



dimension, averaged over many different snapshots. The 
averages are taken with the center of mass of the snap- 
shots shifted to the middle of the simulation box, with 
the constraint that the nematic phase is also located in 
the middle, in accord with Fig. [3 The constraint is nec- 
essary to remove ambiguity arising from cases where the 
isotropic phase is in the middle. Having shifted the cen- 
ter of mass, the density profile is obtained by binning the 
z-axis in steps of Az w 0.17L. The local density p(z) in 
a single snapshot is given by h/vb, with n the number of 
rods in the bin centered around z, and vb the volume of 
a single bin. The density profiles are then averaged over 
all snapshots. Following Ref. Il4t for the bin centered 
around z in a single snapshot, we also define a local ori- 
entational tensor Q(z), calculated using Eq.@ consider- 
ing only the rods inside the bin. The local orientational 
tensor elements are then averaged over all snapshots and 
S(z) — maxev (Q(z)). 

The averaged profiles are shown in Fig. |3| The 
solid curves are hyperbolic tangent fits of the form 
A + Stanh ( z ^ ) Zc ), which describe the data well. Note 
that the profiles are shifted with respect to each other. 
The magnitude of the shift, measured between the inflec- 
tion points, equals S — 0.37 ± 0.04 L. This is consistent 
with theoretical predictions S = 0.45 — 0.5 L 0, as 
well as S pa 0.33 L obtained in simulations of ellipsoids 
|36|. Note that the simulated profiles are broadened due 
to capillary waves [3j. Moreover, we observed consider- 
able fluctuations in the amount of isotropic and nematic 
phase during the simulation, leading to large fluctuations 
in the interface positions along the elongated L z dimen- 
sion. The width of the averaged profile obtained by fix- 
ing the center of mass is therefore additionally broadened 
[37II38I ]. Because of these effects, we cannot compare the 
interfacial width of the simulated profiles to theoretical 
predictions. More important for our purposes, however, 
is the observation that the interfaces are stable. For hard 
rods, the current system size thus seems sufficient to ac- 
commodate stable interfaces. 
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FIG. 3: Nematic order parameter profile S(z) (top), and 
density profile (bottom) for hard rods with L/D — 15 across 
the IN interface. The shift between the profiles is marked S. 
Points are raw simulation data; curves are hyperbolic tangent 
fits. The horizontal lines in the lower frame represent the 
bulk isotropic and nematic densities obtained in the grand 
canonical simulations of Sec. IIV El 



B. Comparison of p and S biased sampling 
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FIG. 4: Coexistence distribution P(p*) for soft rods with 
L/D = 15 and /3e = 2 using box dimensions L x = L y =2.1 L 
and L z = 8.4 L. The solid curve was obtained using a bias on 
the nematic order parameter S; the dashed curve by using a 
bias on the density. The average peak height Af2, multiplied 
by ksT, equals the free energy barrier separating the isotropic 
(ISO) from the nematic (NEM) phase. 



Having established the typical system size required to 
observe stable interfaces, biased sampling on the nematic 
order parameter is considered next. First, we show that 
density and nematic order biased sampling yield the same 
distribution P(N). To this end, we consider a small sys- 
tem of soft rods with L/D — 15 and (3e = 2, in a simula- 
tion box of size L x = L y = 2.1 L and L z = 8.4 L. The lat- 
ter system was investigated in previous work using den- 
sity biased sampling 4j. The corresponding coexistence 
chemical potential reads as /3/i « 5.15. The nematic or- 
der biased sampling scheme is applied to the same system 
using the latter chemical potential and AS = 0.002, see 
Fig. 01 Shown is the coexistence distribution P{N) ob- 
tained using a bias on S (solid curve), as well as using a 
bias on density (dashed curve, reproduced from Ref. |4j). 
The agreement between both methods is strikingly con- 
firmed, thereby justifying the approach of Sec. lIIIBl For 
small systems, the required CPU time is roughly equal 
for both methods. The data sets of Fig.0]required ca. 700 
CPU hours each, on 2.2 GHz Pentium machines. 



C. Interfacial tension of soft rods 

Next, we consider soft rods with L/D = 10 and f3e = 2. 
We aim to accurately measure the interfacial tension. To 
this end, large system sizes are required such that a bias 
on S is essential. As explained before, the elongated L z 
dimension of the simulation box must be large enough 



to accommodate non-interacting interfaces. At the same 
time, L x and L y must be large enough to suppress finite 
size effects in the lateral dimensions. We therefore con- 
sider two system sizes: L x = L y = 3.5 L; L z = 10.5 L 
(system A), and L K — L y = 4L; L z = 14 L (system 
B), where the lateral dimensions are deliberately chosen 
to exceed those of Sec. IIV Al The simulations are per- 
formed using AS = 0.001 and 0.002, for system A and 
B, respectively. An initial estimate of the coexistence 
chemical potential was taken from previous work Q. 

In Fig. [3 the dependence of the nematic order param- 
eter on the number of particle is shown, calculated using 

S(N) = C J2 SiPi(N)e Wi , (8) 

»=1 

with Si = iAS/2 — AS /4, and the remaining symbols 
defined as before. Analogous to fluid-vapor transitions 
[l8L l39l |40| , five distinct regions can be distinguished. In 
region I, a single isotropic phase is observed. Region II 
corresponds to the transition from the bulk isotropic 
phase, to the phase with two parallel interfaces. The 
transition is characterized by the formation of a nematic 
droplet in an isotropic background, which grows with the 
density until it self-interacts through the periodic bound- 
aries, ultimately leading to two parallel interfaces. In re- 
gion III, the interfaces have formed and the system is at 
coexistence, schematically resembling Fig. |5J Increasing 
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FIG. 5: Dependence of the nematic order parameter on the 
density for soft rods with L/D — 10 and fie = 2, obtained 
using two different system sizes. 
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FIG. 7: Logarithm of P(p*) at coexistence for soft rods with 
L/D — 10 and /3e — 2 using two different system sizes. Note 
the fiat region in between the peaks. The arrow indicates 
transition II of Fig. [£] 
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FIG. 6: Weight function W(S) obtained by biasing on the 
nematic order parameter S (see Sec. MI Bl for soft rods with 
L/D = 10 and j3e — 2 using two different system sizes. 



the density further leads to a growth of the nematic do- 
main, at the expense of the isotropic domain. Region IV 
corresponds to the transition to the pure nematic phase, 
during which the system is characterized by an isotropic 
droplet in a nematic background. In region V, finally, a 
single nematic phase is observed. 

In Fig. we show the corresponding weight function 
W(S) for both systems. The double-peaked structure is 



clearly visible. Note that the isotropic peak is signifi- 
cantly higher than the nematic peak. This indicates that 
the chemical potential used in the simulations is below 
the coexistence value. Since coexistence is defined by 
equal weight in the peaks of P(N), and not in W(S), 
Fig. H3 cannot be used to obtain the coexistence chemi- 
cal potential. Instead, P(N) must be constructed first, 
by combining W(S) with the single window distributions 
pi(N); Eq.Q may then be used to extrapolate P(N) to 
coexistence. The resulting coexistence chemical poten- 
tial equals ~ 7.13 for both systems. In Fig. the 
logarithm of P(p*) at coexistence is plotted for both sys- 
tems, scaled with L 2 /(2L 2 ), and the plateaus shifted to 
zero. In this way, the barrier directly reflects the interfa- 
cial tension 7in, in units of k B T/L 2 23] . An important 
observation is that the peaks in both distributions are 
separated by a pronounced flat region. This shows that 
the elongated L z dimension of the simulation box is suffi- 
cient. Moreover, the peak heights are similar, indicating 
that finite size effects in the lateral dimensions L x and L y 
are also small. Therefore, we conclude that the barrier 
in Fig. accurately reflects the interfacial tension 71N 
for soft rods with L/D = 10 and j3e = 2. The resulting 
estimate reads as 7m = 0A9k B T/L 2 = 0.0049 k B T/D 2 . 

In the nematic phase, ca. 6000 rods were simulated for 
system A, and 10,000 for system B. To obtain reliable 
results, a substantial investment in CPU time is thus re- 
quired (ca. 3200 CPU hours were invested for system B). 
Since biased sampling schemes are easy to parallelize, 
results can typically be obtained within 1-2 weeks on a 
modern computer cluster. Accurate sampling is espe- 
cially important around transitions II and IV, and this 
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FIG. 8: Interfacial tension ym(Lx) obtained in cubic systems 
with edge L x , as function of L x , for soft rods with L/D — 10 
and /3e = 2. 



becomes increasingly difficult in large systems 40] . This 
may already be inferred from the scatter in the data of 
system B around transition II (arrow in Fig. 01 • Transi- 
tion IV, on the other hand, is sampled with surprisingly 
little difficulty. The likely explanation is that process II 
requires the formation of a nematic nucleus whose direc- 
tor is aligned in the xy-plane. Process IV, on the other 
hand, does not require any preferred orientation of the 
(isotropic) nucleus, and is therefore easier to sample. 



D. Consequences for finite size extrapolation 

An alternative method to obtain the interfacial tension 
is to measure tin(^x) in cubic systems with edge L x , and 
use the extrapolation equation of Binder |2^| 



7in(£ x ) = Tin + a/L>l + bln(L x )/Ll 



(9) 



to estimate tin. In principle, this approach enables es- 
timates of 7in through an elimination of finite size ef- 
fects, but it requires estimates over a range of values for 
which 7iN(i x )i x /fcBT' S> 1- In practice, however, one 
often tries to use Eq.© using data from smaller sys- 
tems. In Ref. 0, this approach was applied to soft rods 
with L/D = 10 and (3e = 2, assuming b = in Eq.@, 
leading to 7i N = 0.0035 k B T/D 2 . This estimate differs 
profoundly from the one of the previous section, implying 
that finite size extrapolation must be used with care. The 
issue is investigated further in Fig. El Shown is tin(-^x) 
as function of L x , where the open squares are data from 
Ref. 0, and closed squares data from larger systems ob- 
tained in this work. The horizontal line corresponds to 
the estimate of Fig. [7| Note that the data indeed ap- 
proach the latter estimate. The curve is a fit to the open 
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FIG. 9: Coexistence distribution (solid curve) for hard rods 
with L/D = 15, obtained using box dimensions L x = L y = 
10L/3 and L z = 10 L. The barrier Ail is measured between 
the horizontal lines, where the bar gives an indication of the 
uncertainty. The dashed curve shows P(p*) obtained in a 
smaller cubic box with L x = 2.3 L. 



squares using Eq.© with 6 = 0, which summarizes the 
result of Ref. 0. Clearly, the fit fails to capture the data 
of the larger systems. Allowing b in Eq.@ to be non-zero 
will obviously lead to a better fit, but the resulting 7in 
depends sensitively on the range over which the fit is per- 
formed, making this approach somewhat arbitrary. The 
problem partly stems from the difficulty in distinguish- 
ing a/L\ numerically from b ln(L x )/L x , since the range in 
L x that can be sampled is rather small. Additionally, in 
small systems, the interface interactions may be strong. 
This will introduce corrections to Eq.@, which may even 

yield non-monotonic behavior in 71N (£ x ) mmm. a s 

a result, it is difficult to extract Tin vi a finite size extrap- 
olation. In contrast, by using an elongated simulation 
box, and by explicitly checking for finite size effects in 
the lateral dimensions, 71N can be extracted reliably as 
shown in Fig. [7] This, we conclude in hindsight, should 
be the method of choice. 



E. Interfacial tension of hard rods 

Finally, we apply nematic order biased sampling to 
a system of hard rods with L/D = 15, system size 
L x = L y = 10L/3 and L z = 10 L, corresponding to 
ca. 6000 rods in the nematic phase. An initial esti- 
mate of the coexistence chemical potential was obtained 
via Widom insertion |3lj . The nematic order parame- 
ter is sampled with resolution AS = 0.0025 to obtain 
W(S). Combining W(S) with the single window distri- 
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butions pi(N) and applying Eq.JTJ yields for the coexis- 
tence chemical potential /3/i ps 5.58. The corresponding 
coexistence distribution is shown in Fig. El 

Note that P(N) for hard rods is prone to substantial 
statistical error. This is to be expected because the ac- 
ceptance rate of grand canonical insertion for hard rods 
is only 0.004%, compared to 8% for soft rods. Neverthe- 
less, the double-peaked structure is clearly visible. From 
the average peak locations, we obtain pj SO = 0.193 and 
Pnem = 0.220. The latter densities are consistent with 
the bulk plateaus in the density profile, indicated by the 
horizontal lines in Fig. [21 To further check the consis- 
tency of our results, an additional simulation in a smaller 
cubic system with L x = 2.3 L was performed; the cor- 
responding coexistence distribution is shown dashed in 
Fig. El Of course, this system is too small to extract 
the interfacial tension, but the peak positions, and hence 
the coexistence densities, agree well with those of the 
larger system. The agreement with bulk densities ob- 
tained via Gibbs ensemble simulations 0] and Gibbs- 
Duhcm integration [l3| is better than 4%. The height 
of the free energy barrier of the larger system reads 
as Aft = 32 ± 3 k B T, leading to an interfacial tension 
7in ~ 1.4 k B T/L 2 = 0.0064 k B T/D 2 = 0.096 k B T/LD = 
0.10 k B T/(L + D)D. 



V. DISCUSSION AND SUMMARY 

In this paper, we have presented methodic develop- 
ments that allow for the estimation of the interfacial 
tension between isotropic and nematic phases in suspen- 
sions of rods. The problem is challenging because tin 
is very small, and methods that work well for interfaces 
between isotropic phases become problematic, such as 
exploiting the anisotropy of the pressure tensor |33l | . or 
analyzing the capillary wave spectrum (the latter requires 
very precise data from huge systems Q). The novelty of 
the present approach is to combine grand canonical MC 
simulations with a bias on the nematic order parame- 
ter, and obtain tin from the grand canonical distribution 
P(N). The advantage is that the problem of "jamming" 
is largely solved, enabling simulations of large systems. 

The current approach also allows for grand canonical 
simulations of hard rods, enabling a direct comparison 
to theory. In the Onsager limit of infinite rod length, 
theoretical estimates of tin typically range from 0.156 
[TH to 0.34 0, in units of k B T/LD. As expected, this 
exceeds the value for hard rods obtained in this work 
(Tin ~ 0.096 k B T/LD) because L/D = 15 is still far 
from the Onsager limit. As shown by experiment [l| and 
theory [9j, tin increases with L/D. The latter the ory 
is based on the Somoza-Tarazona density functional |45| 
and its main findings are summarized in Fig. 1101 Shown 
are the coexistence densities (left axis) and the interfa- 
cial tension (right axis) as function of the rod elongation 
L/D, where we have adopted the units of Ref. U Open 
and closed squares show the theoretical density of the 




Q 
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FIG. 10: IN coexistence properties of the hard-rod fluid ob- 
tained by theory |j| , compared to simulation results obtained 
in this work. Symbols are explained in the text. 



isotropic and the nematic phase, respectively; the star 
and the cross are the corresponding simulation estimates 
of this work. Closed triangles are the theoretical interfa- 
cial tension, where the line serves to guide the eye; the 
open triangle represents the simulation estimate of tin- 
Theoretical estimates are reported for L/D = 5,10,20, 
but unfortunately not for L/D = 15. This makes a di- 
rect comparison difficult; interpolation of the theoretical 
results, however, seems in good agreement with our sim- 
ulation results, as may be inferred from Fig. 

A typical rod dimension in experiments is L = 150 nm 
and tin = 0.00083 mN/m Q. For T = 298 K, this length 
translates into 0.00025 mN/m using our estimate of Tin- 
Obviously, this estimate differs from the experimental 
one because the hard-rod fluid is a simplified model, but 
it is reassuring to see that the order of magnitude is con- 
firmed. 

The current biased sampling scheme thus seems well 
suited to simulate IN coexistence, even for hard interac- 
tions. Our scheme may also be useful for the application 
of transition path sampling |46| to anisotropic colloidal 
systems, since it can provide valuable starting paths; 
work along these lines is in progress. The remaining 
bottleneck is the low acceptance rate of grand canoni- 
cal insertion. It remains a challenge to address this final 
problem. Since the overall density around the IN transi- 
tion is low, it is anticipated that higher acceptance rates 
can be realized using smarter insertion schemes. To de- 
velop such schemes would be the subject of future work. 
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